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Abstract 


In this article, singularly perturbed differential difference equations having 
delay and advance in the reaction terms are considered. The highest-order 
derivative term of the equation is multiplied by a perturbation parameter 
e taking arbitrary values in the interval (0, 1]. For the small value of ¢, the 
solution of the equation exhibits a boundary layer on the left or right side 
of the domain depending on the sign of the convective term. The terms 
with the shifts are approximated by using the Taylor series approximation. 
The resulting singularly perturbed boundary value problem is solved using 
an exponentially fitted tension spline method. The stability and uniform 
convergence of the scheme are discussed and proved. Numerical exam- 
ples are considered for validating the theoretical analysis of the scheme. 
The developed scheme gives an accurate result with linear order uniform 
convergence. 
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1 Introduction 


A large number of mathematical models have appeared in different areas of 
science and engineering that take into account not just the present state of 
a physical system but also its past history. These models are described by 
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certain classes of functional differential equations often called delay differen- 
tial equations or differential difference equations (DDEs). The class of DDEs 
with characteristics of delay/advance and singularly perturbed behavior is 
known as singularly perturbed differential difference equations (SPDDEs). 
The DDEs with delay or advance term play an important role in model- 
ing many real life phenomena in bioscience, control theory, economics, and 
engineering [5]. Some applications are the mathematical modeling of popu- 
lation dynamics and epidemiology [19], physiological kinetics [4], blood cell 
production [25], and so on. 

In SPDDEs model process, the evaluation not only depends on the current 
state of the system but also includes the past history. A number of model 
problems in science and engineering take the forms of SPDDEs [32]; we list a 
few of them: neuron variability model in computational neuroscience, optimal 
control theory problems, and model describing the motion of sunflower. 

For the perturbation parameter, when € tends to zero, the smoothness of 
the solution of the singularly perturbed problems deteriorates and it forms 
boundary layer; see [6, 27]. In the case where ¢ is very small, standard numer- 
ical methods such as FDM, FEM, and collocation method lead to oscillations 
in the computed solutions. To handle the oscillation, a large number of mesh 
points are required, which is not practical; see [32]. 

The solution methods of SPDDEs have received great attention in recent 
years because of their wide applications. It is of theoretical and practical 
interest to consider numerical methods for such problems [40]. Adilaxmi 
et al. [1, 2] proposed the exponentially fitted nonstandard FDM and in- 
tegration method using a nonpolynomial interpolating function. In articles 
(21, 22, 23, 24], Lange and Miura developed asymptotic methods for solving 
a class of SPDDEs. The authors extend the matched asymptotic method 
initially developed for solving BVPs to obtain an approximate solution for 
SPDDEs. In articles [9, 11, 12, 13], Kadalbajoo and Sharma developed e- 
uniform numerical methods using fitted mesh techniques. Swamy, Phaneen- 
dra, and Reddy [35] used the exponentially fitted Galerkin method for treat- 
ing the problem. The authors in [36, 37] developed a fourth-order FDM with 
an exponential fitting factor. Melesse, Tiruneh, and Deresecite [26] used the 
initial value technique to treat the problem. They showed the applicability 
of the scheme by considering different examples. Ranjan and Prasad [31] 
used the modified fitted FDM for solving the problem. Sirisha, Phaneendra, 
and Reddy [34] developed a finite difference scheme using the procedure of 
domain decomposition. Kumar and Sharma, [20] applied the B-spline collo- 
cation method to approximate the solution of the SPDDEs. Mohapatra and 
Natesan [28] applied the fitted mesh FDM using the equidistributed grid tech- 
nique. In [40], Woldaregay and Duressa developed the exponentially fitted 
FDM with Richardson extrapolation techniques. 

Different authors in [3, 7, 8, 14, 15, 16, 17, 30] have applied the ten- 
sion spline method for treating singularly perturbed reaction diffusion or 
convection diffusion problems. To the best of the authors’ knowledge, the 
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exponentially fitted tension spline method has not been developed for treat- 
ing SPDDEs. Developing uniformly convergent schemes is an active research 
area [33]. This motivates us to develop an accurate and uniformly convergent 
scheme using the exponentially fitted tension spline method. 


Notations: N denotes the number of mesh interval in the discretization, 
C' denotes a positive constant independent of ¢ and N, and the norm ||-|| 
denotes the supremum norm. 


2 Continuous problem 


Consider a class of SPDDE of the form 
—eu" (x) +a(x)u' (x) +a(x)u(x—6)+u(r)u(x)+B8(x)u(x+n) = f(x), «EQ, 


with the interval conditions 


u(x) =9(x), x € [-4,0], 
u(x) =(x), LE el eal, 


where 2. = (0,1), € € (0,1] is the singular perturbation parameter, and 6 
and 7 are delay and advance parameters satisfying 6,7 < ¢. The functions 
a(x), a(x), w(x), B(x), f(z), d(x), and w(x) are assumed to be sufficiently 
smooth and bounded for the existence of unique solution. The coefficient 
functions a(x),w(a), and 8(a) are assumed to satisfy 


a(x) + w(x) + B(2)>atw+B=:0>0, forallzeQ, 


where the constants a,w, and 6 are lower bounds of a(x),w(x), and f(x), 
respectively. 


In the case when 6,7 = 0, equations (1)—(2) reduce to a singularly per- 
turbed boundary value problem, in which for small ¢, it exhibits boundary 
layer. The layer is maintained for 46,7 4 0 but sufficiently small. 


2.1 Properties of the continuous solution 


In the case when 6,7 < €, using Taylor’s series approximation for the terms 
with deviating, the argument is appropriate [38]. Using the Taylor series 
approximation, we approximate 


u(x — 6) u(x) — bu! (x) + (67/2)u" (x) + O(83), 
u(x +) u(x) + nu'(x) + (1°/2)ul(x) + O(n°). 


(3) 
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Replacing (3) in (1) gives 
— ceu" (x) + p(a)u'(x) + g(z)u(x) = f(z), rEQ, (4) 


with the boundary conditions 
u(0) = 6(0), wl) = ¥(1), (5) 


where c. = €? — (3?/2)a— (1?/2)8, p(x) = a(x) — da(a) +6 (x), and q(x) = 
a(x) + B(a) + w(x). For small values of 6 and n, (1)-(2) and (4)-(5) are 
asymptotically equivalent, since the difference between these two equations is 
O(63, 7°). The differential operator L is denoted for the differential equation 
in (4) and defined as 


Lu(x) = —ceul (x) + p(a)u' (x) + q(x)u(z). 


The problem in (4)—(5) exhibits the regular boundary layer of thickness 
O(c-), and the position of the boundary layer depends on the conditions: If 
p(x) <0, then the left boundary layer exists, and if p(a) > 0, then the right 
boundary layer exists. In the case when p(x), x € Q, the change sign interior 
layer will exist [10]. 

The problem obtained by setting c; = 0 in (4)-(5) is called the reduced 
problem and given as 


p(x)ug(x) + q(x)uo(x) = f(x), foralleeQ, 
uo(0) =0(0), uo(1) (1). 


For the right boundary layer case, it does not satisfy the right boundary 
condition, and for the left boundary layer case, it is given as 


(6) 


p(x)ug(x) + g(x)uo(x) = f(x), forallxeQ, 
uo(0) £G(0), wo(1) = W(1). 
For small values of c., the solution u(a) of (4)—(5) is very close to the solution 
uo(x) of (6) or (7). 
Lemma 1 (The maximum principle [40]). For sufficiently smooth function 
z on Q, satisfying z(0) > 0, z(1) = 0, and Lz(xr) > 0, for all x € Q, implies 
that z(a) > 0, for all e € . 


(7) 
( 


Lemma 2 (Stability estimate). The solution u(x) of the continuous equation 
(4)-(5) satisfies the bounded 


|u(x)| < O~* || fll + max{|o(0)|, 'b(1) |}. (8) 


Proof. It is proved by the construction of barrier function and using the 
maximum principle. Let us define barrier functions J+(x) as J=(x) = 
6" || f || + max{(0), 2(1)} + u(x). On the boundary points, we obtain 
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d*(0) = 0 *||f||+ max{o(0), Y(1)} + uO) = 0, 
*(1) = 07? || f+ max{9(0), Y()} + u(1) > 0. 


On the differential operator, we have 


S 


LY* (x) = — ce0'L (x) + p(x)04 (x) + q(x) 04 (2) 

=—c¢.(04u"(x)) + p(x) (04 u'(x)) + q(x) (0 *|| f | + max{(0), b(1)} 
Ze u(x)) 
=q(x) (0 *|| fl + max{9(0), Y(1)}) + f(z) 


>0, since q(x) > 6 > 0. 


By using the hypothesis of the maximum principle, we obtain J*(x) > 
0, for all 2 € Q, which implies the required bound. 


In the next lemma, we obtain a bound for the derivatives of solution. 


Lemma 3. Derivatives of the solutions of the problem in (4)—(5) satisfy the 
bound 


cea rEeQ, 0<k<4, 


Ju*) (a) < C(1 + G* exp ( c 
4 


for the left boundary layer problem and 


rE 


ap ils *)y) 


|u*) (a)| < co(1 a a exp ( 7 
€ 


9 


for the right boundary layer problem. 


Proof. See [6]. 


3 Numerical scheme formulation 


In this article, an exponentially fitted tension spline method is proposed for 
solving equations (1)—(2). The exponential fitting factor is used to hinder the 
influence of the perturbation parameter in the boundary layer region. The 
theory of the asymptotic method is used for developing the exponential fitting 
factor. We consider and treat the left and right boundary layer problems 
separately. 
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3.1 Exponentially fitted tension spline method 


Let 0 = 4% < a1 <a <-+-+-< ay =1 bea uniform partition for [0,1] such 
that 2; = ih, i=0,1,2,...,N. A function S(2,7) = S(z) isa class of C?(Q), 
which interpolates u(x) at the mesh points x; depending on the parameter 
T, reduces to cubic spline in Q as t - 0, and is termed as a parametric cubic 
spline function [17, 3]. In [x;,2;4,], the spline function S(x) satisfies the 
differential equation 


+ [8"(i41) — 79(ti4s) 


(9) 
where S(x;) = u(a;) and tT > 0 is termed as a cubic spline in compres- 
sion. Solving the linear second-order differential equation in (9) and deter- 
mining the arbitrary constants from the interpolation conditions S(aj;.,) = 
u(xi41), 5 (xi) = u(xi), we get 


jt+1 — 2x 


S" (x2) — 7S(x) = [S"(a;) — 7S(2j)] ; 


oo ace. _ [Mita sinh AP) + Mi sinn( A=) y 
h? d? OR pny 2 eee 
-_ xa (Mitt = pz ulwi+i))( a ) T ( 7 AD u(a))(—>=— )], 


where \ = hr'/? and M; =u" (a;) for j =i+1,1. 


Now, differentiating (10) and letting x > 2;, on the interval [x;, 741], we 
obtain 


r 
i na i hA=1)ls 
Mi41(1 Spa? 1 Mécoth A )] 
(11) 


and on the interval [x;_1, x;], we obtain 


S'(az) _ u(x) aa i a [Mi(Acoth A — 1) + Mi-1(1- anh 


Equating the left- and right-hand derivatives at x; gives 


u(x;)—u(ti-1) hy, . A 
; y2 LMi(Acoth A— 1) + Mia(l— =) 


_ mre) u(a;) “ + M;(AcothA — 1)], 


(13) 


r 
Mj.1(1 —- ——~ 
+1 ey 
i=1,2,...,N—1. 


Rearranging, we obtain the tridiagonal system 
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u(xi-1) = 2u(2;) + u(@i+1) 
h2 


AM; -14+202Mi +A, Mia = ; = Me? 


20s IN, 
(14) 
where Ay = x¢(=4 5 — 1) and Ag = ¥4(1 — Acoth)). 
The condition of continuity given in (14) ensures the continuity of the first- 
order derivatives of the spline S'(x) at interior nodes. 
Now, substituting —c-M,; = f(x;) — p(a,;)u'(#;) — q(z;)u(x;) for 7 =1- 1,4 
and 7+ 1, we obtain 
c 
Lu(a;) =— oll 
+ q(i-1)u(ri-1)] + 2A2[p(ai)u' (xi) + o(xi)u(za)] 
+ 1 [p(wig)u' (arias) + q(xi41)U(2i41)| 
=Ai f (@i-1) + 2A2f (wi) + Af (@i41) + Ti(h), t= 1,2,...,N—-1, 
(15) 


u(xj—1) — 2u(xs) + u(ai41)] + Ar[p(ai-1)u’ (w-1) 


where T\(h) is the truncation error in the above discretization, one can see 
the detail in [3], and it is given by 


ht ie 
Ti(h) = (2a + d2)p(wi)u" (Gi) + pt- 12A1)p(wi)ceu (Gi) + O(R*), 
(16) 
for any choice of A; and Ay whose sum is 1/2, except Ay = 1/12 and Ay = 5/12. 
For the choice Ay = 1/12, Az = 5/12, we have 


c-h® (6) 
Ti(h) = 240 & (Gi), Ce © [wi-1, i421]. (17) 


Next, we use the left-shifted, central, and right-shifted finite difference ap- 
proximation as 


—3u(a;-1) + 4u(a;) — u(vi4i) 


Vie. = 
uw (@;-1) Dh + O(h), 
u (ai) Sie +O(h?), and (18) 
2g) ee ee) 008, 
2h 
Substituting (18) in (15) leads to 
e —3u;_-1 + 4u; -— uj; 
Lu; = 7 [usa — 2ug + Us4i] + Ar[p(zi-1)( a“ on +t) + 9(x;_1)uj_1] 
Uj41 — Ui— 3ujs1 — 4uy + Uj 
+ 2A2[p(2i)(>_—=) + q(wi)ui] + Ar [p(@i41)(—= ~) 


2h Qh 
+ q(xig1)uita] = Arf(vi-1) + 2A2f (i) + Af (i421) + To(h), 
a ee eee? 


(19) 
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where u; denotes the approximation of u(x;) in the above discretization and 
To(h) = O(h) + T,(h). 


3.1.1 Case I: Left boundary layer problem 


In this case, the boundary layer occurs on the left side of the domain. From 
the theory of singular perturbations, the zeros-order asymptotic solution of 
(4)—(5) is given as [29] 


u(x) = rig (2) +o”) (4(0)—u9(0)) exp ( [ ( P(z) 12) an) tO(ce). (20) 


ce(x) p(x) 
Using the Taylor series about x = 0 for p(x) and q(x) and simplifying give 
u(x) = u(x) + (4(0) — uo(0)) exp(—p(0)x) + O(ce), (21) 


where uo is the solution of the reduced problem. The domain [0,1] is dis- 
cretized into N equal number of subintervals, each of length h. Let 0 = 29 < 
Ly <@2<+++<a“y =1 be the points such that x; = ih, 1 =0,1,2,...,N. 


Considering h small enough, the discretized form of (21) becomes 
u(ih) ~ uj = uo(th) + ($(0) — wo(0)) exp(—p(0) (ép)), (22) 

where p = h/c- and h =1/N. Similarly, we write 
Uit1 =Uo((é + 1)h) + ($(0) — uo(0)) exp(—p(0)((é + 1)p)), 
Ui-1 =Uo((i — 1)h) + ($(0) — wo(0)) exp(—p(0) ((é — 1)p)). 


In order to handle the influence of the perturbation parameter, the exponen- 
tially fitting factor a, is multiplied on the term containing c, as 


(23) 


CeO —suj—1 + 4uj — Wis 


Du, = h2 [wi-a — 2; + Ui4a] + Ar [p(@i-1)( Dh ) 
+ g(ai—1)t—1] + 2Aalp(ai)( AS) + (waa) 
3uj+4, — 4u; + Uj 
+ Ai [p(wi41)(—> ah =) + q(aiga uit] 
=i f (@i-1) + 2A2f (ai) + Ar f(wi41) + T3(h), 2=1,2,...,N—1. 
(24) 


Multiplying both sides of (24) by h, denoting c:/h = p, and taking the limit 
as h > 0, we obtain 
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p h0 


h-0 
+ 2X2 lim [p(th) (ui4a = Ui-1)| + At lim [p((i + 1)h)(3uj41 _ 4uj; + Ui-1)| = 0. 
h-0 h-0 
(25) 
The expression in (25) simplifies to 
<<, Set 
lim <u — 2uj + isi] = jim (1 + A2)|[p(O)(ui41 — usi-1)]- (26) 
Using the results in (22) and (23), we obtain 
lim [ui+1 — 2t4i + wi-1] =(6(0) — wo(0)) exp(—p(0)ip) 


x [exp(—p(0)p) — 2 + exp(p(0)p)], 
[wit+1 — Ui-1] =((0) — uo(0)) exp(—p(0)ip) 


x [exp(—p(0)p) — exp(p(0)p)]. 


(27) 


lim 

h—0 

Using the result in (27) into (26), we obtain the exponential fitting factor as 
01 = p(0)p(1 + 2) coth(p(0)$). (28) 


2 


Hence, the required finite difference scheme becomes 


CeO —3u;_ + Au; — Uji 
Du, = - = > [usa — Zee; + Us4a] + Ar[p(wi-1)( * +) 
h 2h 
Uj — Ui-~ 
+ g(aj—1)ui—1] + 22[p(ai)(A> >) + a(ai)ui] 


2h 


3uj41 — 4uz + Uj 
+ A [p(ti41)(—* Dh *) g(a aa] 


=Aif (@i_1) + 2raf (i) + A f(@in1) + Ta(h), 1=1,2,...,N—1, 
(29) 


with the boundary conditions up = ¢(0) and uy = w(1). The bound of 
truncation error T3(h) is derived in Theorem 1. 


3.1.2 Case II: Right boundary layer problem 


In this case, the boundary layer is on the right side of the domain. From 
the theory of singular perturbations, the zeros-order asymptotic solution of 
(4)-(5) is given as [29] 
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Using the Taylor series about « = 1 for p(x) and q(x) and simplifying, we 
obtain 


u(a) = u(x) + (C1) — uo(1)) exp (— p(1)(l1—2)) + O(ce), (31) 


where uo is the solution of the reduced problem. The domain [0,1] is dis- 
cretized into N equal number of subintervals, each of length h. Let 0 = 29 < 
Ly <@2<+++<a“2y =1 be the points such that x; = ih, 7 =0,1,2,...,N. 


Considering h is small enough, the discretized form of (21) becomes 
u(ih) ~ uj = uo(th) + (W(1) — uo(1)) exp(—p(1)(1/ce — tp)), (32) 

where p = h/c- and h =1/N. Similarly, we write 
Uit1 =Uo((i + 1)h) + (WL) — uo(1)) exp(—p(1)(1/ce — (¢ + 1)p)), 
Ui—1 =Uo((é — 1)h) + (W(1) — uo(1)) exp(—p(1)(1/ce — (¢ — 1)p)). 


Using the similar procedure as the left boundary layer case, we obtain the 
exponential fitting factor as 


(33) 


02= p(1)p(Ar + d2) cot PUP), (34) 


Hence, the required finite difference scheme becomes 


—8uj-1 + 4u; — Ui4i 
2h 


) + q(2i)ui) 


Lu; == lui — Qu; + wi4a] + A1[p(ai-1)( 
+ q(ai~1)Ui-1] + 2A2[p(a:)( 


3uja1 — 4uj + Uj 
+ Ar [p(eig1)(— oh *) + q(tis1)tes1] 


=A, f(xi-1) + 2d2f (xi) + Arf (vi41) + T3(h), 4=1,2,...,N—1, 


) 


Uit1 — Ui-1 


2h 


with the boundary conditions up = ¢(0) and uy = ~(1). 


3.2 Stability and uniform convergence 


In this section, we discuss the uniform stability and convergence for the right 
boundary layer problems. Similarly, one can do it for the left boundary layer 
case. First, we prove the discrete comparison principle for the scheme in (35) 
for the existence of the unique discrete solution. 


We observe that the nonzero entries of the coefficient matrix of Lu, are 
given by 
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CeO 3p(xji- Ly Lj 
goa = <s Ax ( pies W) res) + q(ai-1)) -— ea ) 
2¢-09 Qp(ai-1) = 2p(ai41) 
aii = Fa A1( h d= qt), (36) 
CeO Li 3p(x; Ty 
ec ate (Pe 1) na rosa) + q(vi41)) a ae ) 
For each 7 = 1,2,...,N —1 and for arbitrary values of h and c-, we have 


Aii-1 < 0, Ayit1 < 0, and aii > 0. 


Lemma 4 (Discrete comparison principle). Assume that, for a mesh func- 
tion u;, there exists a comparison function v; such that L’u; < L’v,, i = 
1,2,...,N—1 and if ug < vp and uy < vy, then u; < v;, i =0,1,2,...,N. 


mm 


Proof. The matrix associated with operator L’, is of size (N +1) x (N +1) 
and where for i = 1 andi = N —1, the terms involving up and uy have been 
moved to the right-hand side. It is easy to see that the matrix of coefficients 
is diagonally dominant and has nonpositive off-diagonal entries. Hence, the 
matrix is an irreducible M matrix. See the details of proof in [18]. 


Lemma 5 (Discrete uniform stability estimate). The solution of the discrete 
scheme in (29) satisfies the bound 


lui] < O-*||Liuill + max{|uo], wv}. (37) 


Proof. Let r = 6~1||L’.u;|| + max{uo, un}, and define the barrier function 
0; by 0} =r+u;. On the boundary points, we obtain 


a 


VE =rtup =0-'||LRu;|| + max{uo, un} + G(0) > 0, 
0% =rtun =07"||LRu;l| + max{ug, uv} + (1) > 0. 


On the discretized spatial domain x;, 0 <i< N, we obtain 


LROF =\— — [(r + Ui-1) — 2(r x Ui) + (r + Ui+1)] 
+ Ax [p(ai-1)( —38(r + wi-1) + “rs ui) — (r+ tit) 
+ (ai_1)(r = uj_a)] + 2A2[p(2;)( ae war) Ur = wi-1) 


t g(a) (7 £ ua)] + Arlp(wi+1)( ee re ee ) 
+ o(ai4a) (7 = i411) 
=[A1q(wi-1) + 2d2q(ai) + drg(i41)] (0 'Lieusll + max{uo, ux}) 
+ [Ai f(wi-1) + 2A2f (wi) + Mf(wi+1)] 20, since qi > 4 > 0. 


From Lemma 4, we obtain J > 0, for all 2; € ON. Hence, the required 
bound is obtained. 
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Now for z > 0, C, and C2 are constants, and we have 


x z= (hj/Ge)" h? 
< h(z)-1<C d ct 38 
Oe eg ee (z) Oy an hice +1 rs (38) 


giving that 


2 
Ice [p(L)p01 + 2) coth (2? _ 1] D*D-u(zi)| x 2 ||” (ai) , (39) 
2 h+e 
since Aj + Ag < 1/2. We obtain the bound 
" +p- = p(l)p +p- ; 
|c-[u (a;) -oD*D u(x,)]| =|6. [p(1)p(Ai + Ag) coth(—>—) — 1] Dt D-u(z;) 
+c. (u’ (xi) — D* D-u(2i))| 
Ch? 
Jue) + Cech? Jue) 


(40) 


Now, let us denote the right-shifted, central, and left-shifted finite differ- 
ences, respectively, as 


_— Ui-1 _— Au; + 3Us41 


DFu(a;) = ah , D®°u(2;) = a and 
—8ui-1 + 4ui — Ui 
Duce _ 38ui-1 — Ui+1 


Using Taylor’s series approximation, we obtain the bound 
lu! (@i-1) — D’ u(ai-1)|SChllu" (Q) |I, 
|u’ (aj) — D°u(x;)|<Ch? || u’"(¢) ||, and (41) 
|w’ (wig) — D¥u(ai1)|<CAl|u"(¢) 


>) 
where ||w'"(¢)|| = maxzo<ai<ay |u"(xs)| and lu!" (¢)|| = maxao<2:<ey |u” (xs)|- 


The next theorem gives the truncation error bound for the proposed 
scheme. 


Theorem 1. Let u(x;) and u; be the solution of (4)-(5) and (29), respec- 
tively. Then, the following error estimate holds: 


| Lu(z) — Lipus| < Ch(1+ cz? exp ( — a) (42) 


Proof. Consider the truncation error bound in the above discretization 


Beponentialhy Atied teased apling wiethad for sineuleuly pertusbell 4. 273 
| Lu(z) - Lipus| <|Lu(2;) _ Lhu(z,)| -- |b" u(2s) -— Lipus| 
<||Ti(A)|| + |ceu" (ai) — ceo D* D~ u(axi)| 

+ |uj_ — D*u(a;)|+|ui, — D?u(z,)| 


+ lug; — D¥u(ai)|- 


(43) 


Using the bounds in (40), (41), and (43) and using 
Ti(h) < Ch*(—2A1 + A2)||w’” (Gi) || + Ceeh*(1 — 12A1) |[u (GI, 


we obtain 


Ch? 
h+ce 
+ ACh? ||u’”’(¢)|]| + Ch4(—2A1 + Az) ||u’” (Ga) 
+ Ceeh4(1 — 12A1)||u (Gi)|]. 


Ju" (aa) |] + ceCh? |u (@,)|] + ArChllu”(O) | 


Lule) —Leu,| < 
|Lu(a;) — Liu| 


Using the bounds for the derivatives of the solution in Lemma 3 gives 


p* (1~a;) =p* (1~2;) 


(1+ e?el= ce )) +. 0n? [ee (1+ eptel =) 


Ch? 
h 
| Lu(2) — Lpuj| = ie 
=p* (1~«;) 


+ Chr, (1+ cel ce )) + Ch? d2(1+ cz3el 


te) 


) 


+ Ch*(—2d; + A2)(1 + eo =) 


+ Cceh*(1 — 1241)(1+ ete =) 
Ch2 = —p*(1—2;) 7 albad 
<a (it etal Ce )) + OF (ce + oe = )) 


=p* (1~«;) —p*(1-2;) 
Ce ) 


+ Chdx (1+ eel )) + Ch? do(1 + eel ) 
=—p* (1-24) 


+ Ch!(—2A1 + Ao)(1-+ pte) 


—p*(1-2;) 


<A jiledee =} 


=p* (1~2j) 


<Ch(1+ cr3el ee )), since: <c.” Sen": 


Lemma 6. For a fixed number of mesh numbers N and for c. > 0, it holds 


exp (=) _ exp (=20=*) 


lim max =—W—s— lim max 
ce 01<i< N-1 roilte ce01<i<N-1 ce 


=0, (44) 


for m = 1,2,3,..., where x; = ih,h =1/N, for all i=1,2,...,N—1. 


Proof. See [39]. 
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Theorem 2. Let u(x;) and u; be the solution of (4)-(5) and (29), respec- 
tively. Then it satisfies the error bound 


sup ||u(a;) — us| < Ch. (45) 
O0<ce.<1 


Proof. Using the result in Lemma 6 into Theorem 1 and using the result in 
Lemma 5, we obtain the required bound. 


4 Numerical results and discussion 


In this section, we consider examples to illustrate the theoretical analysis of 
the proposed scheme. 


Example 1. Consider the problem from [28] 
—eu" (x) +(1+ el”) ul (x) + xe*u(a — 6) + 2?u(x) + (1—e7*)u(z +n) = —-1 
with interval conditions u(z) = 1, —d <a <0,andu(#) = -1, 1 <a < 147. 
Example 2. Consider the problem from [28] 

—eu" (x) + ul (x) — 2u(a — 6) + 5u(x) — u(x +n) = 0 


with interval conditions u(x) = 1, —d <a < Oand u(x) = -1, 1< #< 147. 
The exact solution is given as 


(1 + em emer = (1 + e™ emer 


u(x) = wi ? 
e 2 —eM1 
where 
— =(-1= 26 +n) + /(—1 = 26 + 0)? — 4(e + 6? + 97? /2) 
2(e + 4? + 9/2) 
_ —(-1- 26 +n) (—1 — 26 +)? — 4(e + 62 + 97/2) 
ve 2(e 4+ 62 + 42/2) 


Example 3. Consider the problem from [28] 
—eu" (x) — u' (x) + 2u(a — 6) + 5u(x) — u(x +) =0 


with interval conditions u(x) = 1, -d6 <a <0 and u(x) =0, 1<a<1+7. 
The exact solution is given as 


eMirtme _ eMitmesT 
u(a) ai em2 — em ’ 


where 


Exponentially fitted tension spline method for singularly perturbed ... 275 


_ =(1+ 26 +9) + (1 + 26 + 0)? — 4(€ — 6? + 92/2) 
Ms 2(e — 52 + 12/2) 
_ —(+26 +7) (1+ 26 + )2 — 4(e — 62 + 2/2) 
- ae — 8 + 12/2) 


Example 4. Consider the problem 


—eu" (x)—(1+exp(—2?))u’ (x)—aru(a—6)—2?u(x)—(1.5—exp(—2))u(z+n) = 1 


with interval conditions u(x) =1, —d <a#<Oand u(x) =1, 1<a<1+7. 


Since, the exact solution of the variable coefficient problems is not known, 
we applied the double mesh technique to calculate the maximum absolute 
error. 

Let U denote the computed solution of the problem on N number of 
mesh points and let U?% denote the computed solution on a double number 
of mesh points 2N by including the mid-points 2;41/2 = soe into the 
mesh points. The maximum absolute error is given by 


No _ N No N 2N 
cae — max|U; = u(xi)|, or em — max|U; — U; |, 
and the ¢-uniform error is calculated using EN = maxz,5,» |EX, ie |. The rate of 
: ‘ : ae N 2N 
convergence of the scheme is calculated using rN; ,, = logy (EN;,,/E23,,), and 


the e-uniform rate of convergence is calculated using rN = log, (E“ /E?%). 


Table 1: Example 1, maximum absolute error of the scheme for \1 = 1/12, 2 = 5/12. 


el N =2° 26 2° 28 29 210 


ou 9.4299e-04 5.6749e-04 3.0865e-04 1.6056e-04 8.1847e-05 4.1316e-05 
7 1.4597e-03 6.2724e-04 =2.9779e-04 = =1.4502e-04 = 7.1636e-05 = 3.5613e-05 
os 1.5547e-03 8.9197e-04 4.0302e-04 1.5441e-04 7.2554e-05 = 3.5885e-05 
D2 2.6082e-03  1.0350e-03 4.1282e-04 2.3469e-04 1.0215e-04  3.9059e-05 
2-8 3.0394e-03  1.5279e-03  6.7558e-04 2.6400e-04 1.0473e-04  5.9369e-05 
2-10 3.0383e-03 1.5588e-03 _7.8897e-04 3.8892e-04 1.7038¢-04 6.6326e-05 
2-12 3.0377e-03 —1.5585e-03 7.8910e-04 3.9701e-04 1.9905e-04 —9.7666e-05 
2-14 3.0376e-03 1.5584e-03 7.8907e-04 3.9699e-04 1.9911e-04 9.9709e-05 
2-16 3.0376e-03 1.5584e-03 7.8906e-04 3.9699e-04 1.9911e-04 9.9708e-05 
2-18 — 3.0376e-03 1.5584e-03 7.8906e-04 3.9699e-04 1.991le-04 9.9708¢e-05 
2-20 3.0376e-03 1.5584e-03  _7.8906e-04 3.9699e-04 1.9911e-04 9.9708e-05 


EN 3.0376e-03  1.5584e-03 7.8906e-04 3.9699e-04 1.9911le-04  9.9708e-05 
rN 0.9629 0.9819 0.9910 0.9955 0.9978 - 


Four examples with their solution exhibiting a boundary layer are con- 
sidered. Examples 1 and 2 have the boundary layer on the right side of the 
domain and Examples 3 and 4 have it on the left side of the domain. For 
the detail, one can observe in Figure 1, the layer formation of the solutions 
for different values of ¢«. In Figure 2, the influence of the delay parameter 
on the solution profile is given by considering different values of the delay 
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Figure 1: Boundary layer formation for different values of ¢ and 6 = 0.6e,7 = 0.5e, (a) 
Example 1, (b) Example 2, (c) Example 3 and (d) Example 4. 
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Figure 2: Solution profile for different values of delay parameter for ¢ = 27~?, (a) 
Example 2, (b) Example 4. 


parameter for ¢ = 2~2. The maximum absolute error of the proposed scheme 
for Ay = 1/12 and Az = 5/12 is given in Tables 1-4, for different values of the 
perturbation parameter ¢. One observes that as « — 0 in each column, the 
maximum absolute error becomes stable and uniform. This indicates that 
the proposed scheme is uniformly convergent. In the last two rows of these 
tables, the uniform error and uniform rate of the convergence of the scheme 
are given. In Table 5, the rate of convergence of the scheme is given for dif- 
ferent values of ¢ ranging from 2~!* to 2~7°. It is observed that the scheme 
gives a linear order uniform convergence. In Tables 6-8, the comparison of 
the proposed scheme with the result in [28] is given. As we observe, the uni- 
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Table 2: Example 2, maximum absolute error of the scheme for Ay = 1/12, A2 = 5/12. 


el N=25 26 a 28 29 gto 


20 7.4098e-04  3.7006e-04 1.8494e-04 9.2446e-05 4.6218e-05 = 2.3108e-05 
ox 2.1452e-03 1.0603e-03 5.2699e-04 2.6277e-04 =1.3121e-04 6.5561e-05 
Q-4 3.3473e-03 1.5615e-03 = 7.5382e-04 = 3.7142e-04 —-1.8438e-04 = 9.1872e-05 
2-6 4.4315e-03  2.1468e-03 9.3461e-04 4.3381le-04 2.0955e-04 1.0317e-04 
a 6.2230e-03  2.7096e-03 =1.1361e-03  5.5377e-04 2.4090e-04 1.1169e-04 
2-10 6.2230e-03 3.2715e-03 -1.6132e-03  6.9087e-04 2.8575e-04 1.3957e-04 
2-12 6.3987e-03 3.2747e-03  1.6565e-03  8.3275e-04 4.0708e-04 1.7359e-04 
2-14 6.4018e-03 3.2751e-03 1.6567e-03 8.3326e-04 4.1787e-04  2.0915e-04 
2-16 6.4025e-03 3.2752e-03  1.6568e-03 8.3329e-04 4.1788e-04 2.0915e-04 
2-18 6.4025e-03 3.2752e-03 1.6568e-03 8.3329e-04 4.1788e-04 2.0915e-04 
2-20 6.4025e-03 3.2752e-03 1.6568e-03 8.3329e-04 4.1788e-04 2.0915e-04 


EN 6.4025e-03  3.2752e-03 = 1.6568e-03 = 8.3329e-04 = 4.1788e-04 2.0915e-04 
rN 0.9671 0.9832 0.9915 0.9957 0.9986 - 


Table 3: Example 3, maximum absolute error of the scheme for Ay = 1/12, A2 = 5/12. 


el N=25 26 Q7 28 29 210 


Q0 9.7173e-04 4.7967e-04 2.3842e-04 1.1886e-04 5.9346e-05 2.9652e-05 
2-2 2.2582e-03  1.0918e-03 5.3781e-04 2.6712e-04 1.331le-04 6.6445e-05 
2-4 4.1684e-03  1.9446e-03 9.3525e-04 4.6062e-04 2.2884e-04 1.1409e-04 
2-8 9.0197e-03  3.4155e-03 = 1.3580e-03 = 3.1052e-04 =. 2.0955e-04_ =—-1.5427e-04 
Q-8 1.3302e-02  6.2951e-03  2.4966e-03 -9.1533e-04 = 3.7245e-04 —-1.7575e-04 
2-10 1.3600e-02 7.2675e-03 3.6950e-03 1.6616e-03 6.4185e-04 2.3339e-04 
2-12 1.3615e-02 —7.2783e-03  3.7679e-03  1.9170e-03 9.4973e-04 4.2126e-04 
2-14 1.3618e-02 7.2804e-03 3.7690e-03  1.9182e-03 9.677le-04 4.8587e-04 
2-16 —1.3619e-02 7.2809e-03 3.7690e-03 -1.9183e-03 9.6778e-04  4.8607e-04 
2-18 — 1.3619e-02 7.2809e-03  3.7690e-03  1.9183e-03 9.6778e-04 4.8607e-04 
2-20 1.3619e-02 7.2809e-03  3.7690e-03 1.9183e-03 9.6778e-04 4.8607e-04 


EN 1.3619e-02 7.2809e-03  3.7690e-03 1.9183e-03 9.6778e-04  4.8607e-04 
rN 0.9034 0.9499 0.9744 0.9871 0.9935 - 


form error and uniform rate of convergence of the proposed scheme is better 
than that of in [28]. 


5 Conclusion 


This article dealt with the numerical treatment of SPDDEs having shifts on 
the reaction terms. The solution of the considered problem exhibited the 
boundary layer on the left or right side of the domain as e« — 0. The terms 
involving the shift were approximated using the Taylor series approximation. 
The exponentially fitted tension spline method was used for treating the 
resulting singularly perturbed boundary value problem. The first derivative 
terms were approximated using left-shifted, central, and right-shifted finite 
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Table 4: Example 4, maximum absolute error of the scheme for Az = 1/12, A2 = 5/12. 


el N=25 26 ar 28 29 gto 

Qo 1.9784e-04 9.7513e-05 4.8442e-05 2.4139e-05 1.2049e-05 —6.0194e-06 
2-2 3.7410e-04 2.0860e-04 =—1.0982e-04 5.6313e-05 2.8509e-05 1.4342e-05 
2-4 1.4674e-04 1.6633e-04 1.2390e-04 7.2396e-05 3.8825e-05  2.0070e-05 
2-6 1.4706e-03 3.3509e-04 4.43897e-05 4.6273e-05 3.5377e-05  2.0784e-05 
2-8 — 2.0794e-03  1.0391e-03  3.9571e-04 8.8140e-05 1.1611e-05 = 1.2172e-05 
2-10 2.0806e-03 1.0780e-03 5.4814e-04 2.6668e-04 1.0074e-04 2.2328e-05 
2-12 2.0806e-03 1.0780e-03 5.4844e-04  2.7658e-04 1.3881e-04 6.7101e-05 
2-14 2.0806e-03 1.0780e-03 5.4844e-04 2.7658e-04 1.3888e-04 6.9589e-05 
2-16  2.0806e-03 1.0780e-03 5.4844e-04 2.7658e-04 1.3888e-04 6.9589e-05 
2-18  2.0806e-03 1.0780e-03 5.4844e-04 2.7658e-04 1.3888e-04 6.9589e-05 
2-20  2.0806e-03 1.0780e-03 5.4844e-04 2.7658e-04 1.3888e-04 6.9589e-05 
EN 2.0806e-03 1.0780e-03 5.4844e-04 2.7658e-04 1.3888e-04 6.9589e-05 
rN 0.9486 0.9750 0.9876 0.9939 0.9969 - 


Table 5: Example 1, Cane 


el N=2° 26 Q7 28 29 
Example 1 
g-12 0.9628 0.9819 0.9910 0.9960 1.0272 
ont 0.9629 0.9819 0.9910 0.9955 0.9978 
2-28 0.9629 0.9819 0.9910 0.9955 0.9978 
Q18 0.9629 0.9819 0.9910 0.9955 0.9978 
Q-20 0.9629 0.9819 0.9910 0.9955 0.9978 
Example 4 
gre 0.9486 0.9750 0.9876 0.9946 1.0487 
o-14 0.9486 0.9750 0.9876 0.9939 0.9969 
o-16 0.9486 0.9750 0.9876 0.9939 0.9969 
g-18 0.9486 0.9750 0.9876 0.9939 0.9969 
g720 0.9486 0.9750 0.9876 0.9939 0.9969 


) of the scheme for Ay = 1/12, A2 = 5/12. 


Table 6: Comparison of uniform error and uniform rate of convergence of Example 1. 


el N =2° 28 oe 28 29 210 
Propose Scheme 

EN — 3.0376e-03  1.5584e-03 7.8906¢-04  3.9699e-04 1.991le-04 9.9708e-05 

rN 0.9629 0.9819 0.9910 0.9955 0.9978 - 
Result in [28] 

EN  8.9743e-02 4.6893e-02 2.4148e-02 1.5602e-02 7.5110e-03 3.5319e-03 

rN 0.9364 0.9575 0.6302 1.0547 1.0886 - 


difference approximation. The stability of the scheme was investigated using 
the comparison principle and solution bound. The uniform convergence of 
the scheme was proved, and it gave the first-order uniform convergent. The 
performance of the scheme was compared with some published articles and 
it gave an accurate result. 
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Table 7: Comparison of uniform error and uniform rate of convergence of Example 2. 


ea). N =2° 26 hs 28 29 210 
Propose Scheme 

EN  6.4025e-03 3.2752e-03 1.6568e-03  8.3329e-04 4.1788e-04 2.0915e-04 

rN 0.9671 0.9832 0.9915 0.9957 0.9986 - 
Result in [28] 

EN  1.0273e-01 6.1537e-02 3.8643e-02 2.2077e-02 1.2395e-02 7.0772e-03 

rN 0.7393 0.6712 0.8074 0.8328 0.8085 - 


Table 8: Comparison of uniform error and uniform rate of convergence of Example 3. 


el N=25 26 a 28 29 210 
Propose Scheme 

EN  1.3619e-02 7.2809e-03 3.7690e-03 1.9183e-03 9.6778e-04 4.8607e-04 

rN 0.9034 0.9499 0.9744 0.9871 0.9935 - 
Result in [28] 

EN  9.6126e-02 5.7165e-02 3.3247e-02 1.8984e-02 1.0685e-02 5.9444e-03 

rN 0.7498 0.7819 0.8084 0.8293 0.8459 - 
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